Teoria Macroeconômica II - Artigo parte II - Sétimo Var

Introdução

Este documento contém a sétima especificação do VAR que eu testei. Ela é a mesma variação do VAR 1, porém com diagnóstico de convergência.

Para as transformações que foram feitas nos dados e a motivação da escolha das variáveis, sugiro ver meu arquivo de descritivas.

Para fazer o BVar utilizei o pacote bvarsv do Fabian Kruger. Em particular, esses dois arquivos me ajudaram muito: bvarsv: An R implementation of the Primiceri (2005) model for macroeconomic time series e Replication of figures in Del Negro and Primiceri (2015).

Para fazer o diagnóstico da cadeia de Markov (algo que não estava presente nos outros VARs), usei o pacote coda junto com …

Download de pacotes e dependências

Se o mirror 10 (UFRJ) estiver fora do ar, tente mudar ind = 10 para ind = 1. Mesma coisa se o BETS não instalar.

chooseCRANmirror(graphics = FALSE, ind = 10)
if (!require("pacman")) install.packages("pacman")
pacman::p_load(ggplot2, forecast, BETS, seasonal, seasonalview, bvarsv, lubridate, zoo, stargazer, gridExtra, reshape2, ggfortify, RColorBrewer, scales, coda, foreach, doParallel, knitr, grid, ggpubr) 

Algumas funções auxiliares:

## Peguei essa função pronta
matplot2 <- function(...){
  matplot(..., type = 'l', lty = 1, lwd = 2, bty = "n", ylab = "")
}

# Essa função calcula os quantis e organiza 1 quantil, média, 3 quantil.
stat.helper <- function(z) c(mean(z), quantile(z, c(0.16, 0.84)))[c(2,1,3)]

# Paleta de cores
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

cols1 <- cbPalette[c(2,4,2)]
cols2 <- cbPalette[c(2,4,6)]

Especificação 1: 5 variáveis, 2003-2017

Tentativa 1

A primeira tentativa foi com dados de janeiro de 2003 a outubro de 2017 e é igual à especificação do primeiro VAR. Minha intenção era testar os diagnósticos e ver se o computador consegue salvar a séries completas do amostrados de Gibbs.

1. Razão capital-trabalho, calculada a partir das seguintes séries:
   • Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (7620)
   • Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (7621)

Passei um filtro para tirar sazonalidade, sem transformação prévia nos dados (ajuste default na função seas exceto no transformation que ficou igual a NONE).

2. Taxa Selic calculada a partir da série:
   • Taxa de juros - Selic acumulada no mês (4390). Passei de taxa mensal para anual usando a fórmula: \(\left((1+tx/100)^12 -1\right)*100\) e passei um filtro para tirar sazonalidade, sem transformação prévia nos dados (ajuste default na função seas exceto no transformation que ficou igual a NONE).
3. IPCA (para inflação), utilizando a seguinte série:
   • Índice nacional de preços ao consumidor-amplo (IPCA) (433). Como ele está em variação percentual mensal, foi necessário transformar para acumulado dos últimos 12 meses utilizando a fórmula: \[IPCA_i = \left[\left(\prod\limits_{j=i-11}^i \left(\frac{IPCA_j}{100}+1\right) \right) -1\right]*100\]
4. IBC-Br (para o PIB per capita), utilizando a seguinte série:
   • Índice de Atividade Econômica do Banco Central (IBC-Br) - com ajuste sazonal (24364). Fiz a diferença da série em log, igual o câmbio.
5. Crescimento da taxa de câmbio (para growth of the nominal effective exchange rate). Usei a série:
   • Taxa de câmbio - Livre - Dólar americano (venda) - Fim de período - mensal (3696). Como o Mumtaz disse que usou log diferença, eu calculei o log e tirei primeira diferença nessa série.

Download dos dados

#rm(list = ls())

# Auxiliary variables, so I don't need to bother when something changes
inicio <- "2003-02-01"
fim    <- "2017-10-31"

# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca   <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries

trabalho <- BETS.get("7620", from = inicio, to = fim)
capital  <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho

# Usando ndiffs(ibcbr,test="adf",alpha = 0.1) se encontra que o ibcbr é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
ibcbr_raw    <- BETS.get("24364", from = inicio_cambio, to = fim)
ibcbr        <- diff(log(ibcbr_raw), 1)

cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio     <- diff(log(cambio_raw), 1)

selic_4390 <- BETS.get("4390", from = inicio, to = fim) 
selic <- ((1+selic_4390/100)^(12)-1)*100

ipca_raw <- BETS.get("433", from = inicio_ipca, to = fim) 
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()

final <- length(ipca_acum)
for (i in 12:final){
  ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}

ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))

ipca <- ts(ipca,  start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD

m <- seas(x = capital_trabalho)
capital_trabalho2 <- final(m)

m <- seas(x = selic, transform.function = "none")
selic2 <- final(m)

df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho2, selic2, ibcbr, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "IBCBr", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")

cores <- brewer.pal(5, "Dark2")

# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[1])+
        scale_y_continuous(name="Capital trabalho") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p2 <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[2])+
        scale_y_continuous(name="Selic (%a.a.)") +
        scale_x_date(date_breaks = "12 months")+ 
        theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p3 <- ggplot(df2[which(df2$variable == "IBCBr"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[3])+
        scale_y_continuous(name="IBC-Br") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[4])+
        scale_y_continuous(name="Tx. Câmbio (var)") +
        scale_x_date(date_breaks = "12 months")+
        theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

p5 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, colour = cores[5])+
        scale_y_continuous(name="IPCA (acum. 12m.)") +
        scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
        theme_bw()
p5 <- p5 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))

grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5)

rm(capital, capital_trabalho, trabalho, ipca_acum, ipca_raw, ibcbr_raw, selic_4390, selic, cambio_raw)
rm(p1, p2, p3, p4, p5, df2)

### Descriptives

descriptives     <- matrix(NA, nrow = 8, ncol = (ncol(df1)-1))
rownames(descriptives) <- c("Observações", "Mínimo", "1o quartil",
                      "Média", "Mediana",  "3o quartil", "Máximo",
                      "Desv. Pad.")

colnames(descriptives) <- names(df1)[-1]

desc <- function(x) {
  n       <- length(x)
  minimum <- min(x, na.rm = TRUE)
  first_q <- quantile(x, 0.25, na.rm = TRUE)
  media   <- mean(x, na.rm = TRUE)
  mediana <- median(x, na.rm = TRUE)
  third_q <- quantile(x, 0.75, na.rm = TRUE)
  maximum <- max(x, na.rm = TRUE)
  std     <- sd(x, na.rm = TRUE)

    return(list(n = n, minimum = minimum, first_quar = first_q, media = media, mediana = mediana, third_quar = third_q, maximum = maximum, std = std))
}

for (i in 1:8){
  descriptives[i, 1] <- round(as.numeric(desc(df1[,2])[i]),4)
  descriptives[i, 2] <- round(as.numeric(desc(df1[,3])[i]),4)
  descriptives[i, 3] <- round(as.numeric(desc(df1[,4])[i]),4)
  descriptives[i, 4] <- round(as.numeric(desc(df1[,5])[i]),4)
  descriptives[i, 5] <- round(as.numeric(desc(df1[,6])[i]),4)
}

descriptives[1,] <- as.integer(descriptives[1,])

descriptives <- data.frame(descriptives)

stargazer(descriptives, summary=FALSE, header = TRUE, type = 'html')

	Capital_trabalho	Selic	IBCBr	Cambio	IPCA
Observações	177	177	177	177	177
Mínimo	0.190	6.707	-0.031	-0.149	2.456
1o quartil	0.436	10.305	-0.003	-0.027	4.832
Média	0.526	12.953	0.002	-0.0004	6.497
Mediana	0.493	12.263	0.002	-0.005	6.037
3o quartil	0.570	14.416	0.007	0.022	7.138
Máximo	1.305	26.608	0.019	0.158	17.235
Desv. Pad.	0.163	4.105	0.008	0.045	2.799

#stargazer(descriptives, summary=FALSE, header = TRUE, type = 'latex')

rm(descriptives, df1)

Ajustando um modelo BVar

Inicialmente ajustei um modelo com uma defasagem e que utiliza as primeiras 24 observações para fazer a estimativa de MQO para jogar na priori.

var1 <- cbind(capital_trabalho2, selic2, ibcbr, cambio, ipca) # The bvar function does not allows data.frames

set.seed(1)

nburn. <- 5000
nrep. <- 50000

fit1 <- bvar.sv.tvp(var1, p=1, tau = 24, nburn = nburn., nrep = nrep., nf = 1, thinfac = 1)

## [1] "2018-01-26 06:55:45 -- now starting MCMC"
## [1] "2018-01-26 06:56:51 -- now at iteration 5000"
## [1] "2018-01-26 06:59:11 -- now at iteration 15000"
## [1] "2018-01-26 07:01:26 -- now at iteration 25000"
## [1] "2018-01-26 07:03:43 -- now at iteration 35000"
## [1] "2018-01-26 07:06:00 -- now at iteration 45000"
## [1] "2018-01-26 07:08:17 -- now at iteration 55000"

Alguns gráficos

## Colocando as datas como strings
tm1 <- as.yearmon(time(var1))

# Eixo x
xax <- time(var1)
# Perde uma observação pelo lag e outras 24 pela amostra para estimar por MQO, então começa no tempo 26
xax <- xax[26:177]

# Marcas verticais
gp <- seq(1999, 2017, 1)

Funções impulso resposta

Peguei a explicação daqui daqui: A FIR estima o impacto de um choque unitário em algum elemento de \(\varepsilon_t\). Os elementos de \(c_t\), \(B_{j,t}\), \(A_t\) e \(\Sigma_t\) em (1) são mantidos fixos em seus valores em \(t\). Note que a ordem das variáveis importa e deve ser justificada por argumentos econômicos. As FIR contém os percentis 5, 25, 50, 75 e 95.

# Função Impulso resposta da Selic na razão capital-trabalho
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 1)

# Função Impulso resposta da Selic no IBC-Br
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 3)

# Função Impulso resposta da Selic no câmbio
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 4)

# Função Impulso resposta Selic no IPCA
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 5)

# Função Impulso resposta do IBC-Br na razão capital trabalho
ira <- impulse.responses(fit1, impulse.variable = 3, response.variable = 1)

# Função Impulso resposta do Cambio na razão capital trabalho
ira <- impulse.responses(fit1, impulse.variable = 4, response.variable = 1)

# Função Impulso resposta do IPCA na razão capital trabalho
ira <- impulse.responses(fit1, impulse.variable = 5, response.variable = 1)

Avaliando os parâmetros

rm(var1, ira)

# Saving parameters
beta_ct     <- parameter.draws(fit1, type = "lag1", row = 1, col = 1)
beta_selic  <- parameter.draws(fit1, type = "lag1", row = 2, col = 2)
beta_ibcbr  <- parameter.draws(fit1, type = "lag1", row = 3, col = 3)
beta_cambio <- parameter.draws(fit1, type = "lag1", row = 4, col = 4)
beta_ipca   <- parameter.draws(fit1, type = "lag1", row = 5, col = 5)

sd_ct       <- parameter.draws(fit1, type = "vcv", row = 1, col = 1)
sd_selic    <- parameter.draws(fit1, type = "vcv", row = 2, col = 2)
sd_ibcbr    <- parameter.draws(fit1, type = "vcv", row = 3, col = 3)
sd_cambio   <- parameter.draws(fit1, type = "vcv", row = 4, col = 4)
sd_ipca     <- parameter.draws(fit1, type = "vcv", row = 5, col = 5)

# salva o fit em um csv e limpa da memória
saveRDS(fit1, "C:\\Users\\aisha\\Documents\\Var da Aisha\\Var7\\fit1.rds")
#fit1 <- readRDS("C:\\Users\\aisha\\Documents\\Var da Aisha\\Var7\\fit1.rds")
rm(fit1)

# Convert stuff into mcmc objects
mcmc_beta_ct       <- mcmc(beta_ct, start = 1, end = 50000)
mcmc_beta_selic    <- mcmc(beta_selic, start = 1, end = 50000)
mcmc_beta_ibcbr    <- mcmc(beta_ibcbr, start = 1, end = 50000)
mcmc_beta_cambio   <- mcmc(beta_cambio, start = 1, end = 50000)
mcmc_beta_ipca     <- mcmc(beta_ipca, start = 1, end = 50000)

mcmc_sd_ct         <- mcmc(sd_ct, start = 1, end = 50000)
mcmc_sd_selic      <- mcmc(sd_selic, start = 1, end = 50000)
mcmc_sd_ibcbr      <- mcmc(sd_ibcbr, start = 1, end = 50000)
mcmc_sd_cambio     <- mcmc(sd_cambio, start = 1, end = 50000)
mcmc_sd_ipca       <- mcmc(sd_ipca, start = 1, end = 50000)

Os gráficos abaixo contém os betas (esq) e volatilidade (dir) ao longo do tempo, com o 16o e o 68o percentis (como foi feito no artigo do Primiceri).

# Gráfico do desvio padrão do resíduo da razão capital trabalho
x1     <- t(apply(sqrt(sd_ct), 2, stat.helper))
x1betact <- t(apply(beta_ct, 2, stat.helper))
#x1beta <- cbind(desc_betact$quantiles[,1], desc_betact$statistics[,1], desc_betact$quantiles[,5]) # ficou muito amplo mas meu coração não deixou eu apagar essa linha
desc_betact <- summary(mcmc_beta_ct)
desc_sdct   <- summary(mcmc_sd_ct)

old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x1betact, ylim = c(min(desc_betact$quantiles[,1]), max(desc_betact$quantiles[,5])), col = cols1, main = "Betas da Razão Capital/Trabalho", xlab = "Tempo")
matplot2(x = xax, y = x1, ylim = c(sqrt(min(desc_sdct$quantiles[,1])), sqrt(max(desc_sdct$quantiles[,5]))), col = cols1, main = "Volatilidade da Razão Capital/Trabalho", xlab = "Tempo")

#matplot2(x = xax, y = x1, ylim = c(0.04, 0.065), col = cols1, main = "Volatilidade da Razão Capital/Trabalho", xlab = "Tempo")
#matplot2(x = xax, y = x1, ylim = c(min(x1[,1])-2*sd(x1[,1]), max(x1[,3])+2*sd(x1[,3])), col = cols1, main = "Volatilidade da Razão Capital/Trabalho", xlab = "Tempo")
par(old.par) # resets par()

rm(desc_betact, desc_sdct)

# Gráfico do desvio padrão do resíduo da selic
x2 <- t(apply(sqrt(sd_selic), 2, stat.helper))
x2betaselic <- t(apply(beta_selic, 2, stat.helper))
desc_betaselic <- summary(mcmc_beta_selic)
desc_sdselic   <- summary(mcmc_sd_selic)

old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x2betaselic, ylim = c(min(desc_betaselic$quantiles[,1]), max(desc_betaselic$quantiles[,5])), col = cols1, main = "Betas da Selic", xlab = "Tempo")
matplot2(x = xax, y = x2, ylim = c(sqrt(min(desc_sdselic$quantiles[,1])), sqrt(max(desc_sdselic$quantiles[,5]))), col = cols1, main = "Volatilidade da Selic", xlab = "Tempo")

par(old.par) # resets par()

rm(desc_betaselic, desc_sdselic)

# Gráfico do desvio padrão do resíduo do ibcbr
x3 <- t(apply(sqrt(sd_ibcbr), 2, stat.helper))
x3betaibcbr <- t(apply(beta_ibcbr, 2, stat.helper))
desc_betaibcbr <- summary(mcmc_beta_ibcbr)
desc_sdibcbr   <- summary(mcmc_sd_ibcbr)

old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x3betaibcbr, ylim = c(min(desc_betaibcbr$quantiles[,1]), max(desc_betaibcbr$quantiles[,5])), col = cols1, main = "Betas do IBC-Br", xlab = "Tempo")
matplot2(x = xax, y = x3, ylim = c(sqrt(min(desc_sdibcbr$quantiles[,1])), sqrt(max(desc_sdibcbr$quantiles[,5]))), col = cols1, main = "Volatilidade do IBC-Br", xlab = "Tempo")

par(old.par) # resets par()

rm(desc_betaibcbr, desc_sdibcbr)

# Gráfico do desvio padrão do resíduo do cambio
x4 <- t(apply(sqrt(sd_cambio), 2, stat.helper))
x4betacambio <- t(apply(beta_cambio, 2, stat.helper))
desc_betacambio <- summary(mcmc_beta_cambio)
desc_sdcambio   <- summary(mcmc_sd_cambio)

old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x4betacambio, ylim = c(min(desc_betacambio$quantiles[,1]), max(desc_betacambio$quantiles[,5])), col = cols1, main = "Betas da Variação da Tx. de Câmbio", xlab = "Tempo")
matplot2(x = xax, y = x4, ylim = c(sqrt(min(desc_sdcambio$quantiles[,1])), sqrt(max(desc_sdcambio$quantiles[,5]))), col = cols1, main = "Volatilidade da Variação da Tx. de Câmbio", xlab = "Tempo")

par(old.par) # resets par()

rm(desc_betacambio, desc_sdcambio)

# Gráfico do desvio padrão do resíduo do ipca
x5 <- t(apply(sqrt(sd_ipca), 2, stat.helper))
x5betaipca <- t(apply(beta_ipca, 2, stat.helper))
desc_betaipca <- summary(mcmc_beta_ipca)
desc_sdipca   <- summary(mcmc_sd_ipca)

old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x5betaipca, ylim = c(min(desc_betaipca$quantiles[,1]), max(desc_betaipca$quantiles[,5])), col = cols1, main = "Betas do IPCA acum. 12 meses", xlab = "Tempo")
matplot2(x = xax, y = x5, ylim = c(sqrt(min(desc_sdipca$quantiles[,1])), sqrt(max(desc_sdipca$quantiles[,5]))), col = cols1, main = "Vol. do IPCA acum. 12 meses", xlab = "Tempo")

par(old.par)

rm(desc_betaipca, desc_sdipca)

Diagnóstico

Estou me baseando nos seguintes linkis: * www.patricklam.org/teaching/convergence_print.pdf para traceplot, density plot, running mean plot (é o plot da iteração contra a média das observações até a iteração (ele está indo de 10 em 10 para frente na média)) * https://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/

Diagnósticos para as médias

Um dos diagnósticos é plotar a cadeia (usando plot(mcmc_object)). Ele retorna o gráfico dos valores (um gráfico de linha - quanto mais ruidoso melhor) e um traceplot que é o gráfico da distribuição. Se eu não estou enganada, os meus Betas tem distribuição à posteriori normal e falta ver a distribuição dos sigmas. Se a densidade das matrizes é inversa wishart eles tem distribuição gamma, mas acho que não tem fórmula fechada para \(\Sigma\) - preciso ver no artigo do Primeri depois.

Outro diagnóstico das médias envolve fazer o gráfico da média usando uma janela móvel para ver se a posteriori cai em alguma região e fica empacada.

Como eu tenho muitos instantes de tempo (e portanto uma quantidade enorme de betas e de sigmas), eu selecionei arbitrariamente 15 betas (que é a divisão inteira de t por 10 - então estou plotando 10% dos parâmetros para cada \(y\)).

sequencia <- seq(10, dim(beta_ct)[2], 10)
sequencia2 <- seq(1:dim(beta_ct)[1])

#####################################
# Para a razão capital trabalho
#####################################

#############
# Betas
#############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_beta_ct[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "da Razão Cap./Trab."))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,mcmc_beta_ct[,i])
  names(df1) <- c("Iteracao", "Beta")
  ggplot(df1, aes(Iteracao, Beta)) +
    geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

##############
# Volatilidade
##############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_sd_ct[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Razão Cap./Trab."))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,sqrt(mcmc_sd_ct[,i]))
  names(df1) <- c("Iteracao", "DP")
  ggplot(df1, aes(Iteracao, DP)) +
    geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

#####################################
# Para a Selic
#####################################

#############
# Betas
#############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_beta_selic[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "da Selic"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,mcmc_beta_selic[,i])
  names(df1) <- c("Iteracao", "Beta")
  ggplot(df1, aes(Iteracao, Beta)) +
    geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

##############
# Volatilidade
##############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_sd_selic[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Selic"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,sqrt(mcmc_sd_selic[,i]))
  names(df1) <- c("Iteracao", "DP")
  ggplot(df1, aes(Iteracao, DP)) +
    geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)


#####################################
# Para o IBC-Br
#####################################

#############
# Betas
#############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_beta_ibcbr[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "do IBC-Br"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,mcmc_beta_ibcbr[,i])
  names(df1) <- c("Iteracao", "Beta")
  ggplot(df1, aes(Iteracao, Beta)) +
    geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

##############
# Volatilidade
##############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_sd_ibcbr[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Selic"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,sqrt(mcmc_sd_ibcbr[,i]))
  names(df1) <- c("Iteracao", "DP")
  ggplot(df1, aes(Iteracao, DP)) +
    geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

#####################################
# Para o Câmbio
#####################################

#############
# Betas
#############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_beta_cambio[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "do Câmbio"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,mcmc_beta_cambio[,i])
  names(df1) <- c("Iteracao", "Beta")
  ggplot(df1, aes(Iteracao, Beta)) +
    geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

##############
# Volatilidade
##############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_sd_cambio[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Selic"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,sqrt(mcmc_sd_cambio[,i]))
  names(df1) <- c("Iteracao", "DP")
  ggplot(df1, aes(Iteracao, DP)) +
    geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)


#####################################
# Para a IPCA
#####################################

#############
# Betas
#############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_beta_ipca[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "do IPCA"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,mcmc_beta_ipca[,i])
  names(df1) <- c("Iteracao", "Beta")
  ggplot(df1, aes(Iteracao, Beta)) +
    geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

##############
# Volatilidade
##############

# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  plot(mcmc_sd_ipca[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "do IPCA"))
}

par(old.par)  

# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  df1 <- data.frame(sequencia2,sqrt(mcmc_sd_ipca[,i]))
  names(df1) <- c("Iteracao", "DP")
  ggplot(df1, aes(Iteracao, DP)) +
    geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
    scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
    theme_bw()
} 
par(old.par)

Agora tem o cumuplot() que faz o gráfico dos quantis da cadeia (de forma acumulada). O default é 0.025, 0.5, 0.975. Se eu tiver tempo depois, vou mexer no layout desses gráficos pois são feiosos demais.

#####################################
# Para a razão capital trabalho
#####################################

#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_beta_ct[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "da Razão Cap./Trab."))
}

#par(old.par) 


#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_sd_ct[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "da Razão Cap./Trab."))
}

#par(old.par) 

#####################################
# Para a selic
#####################################

#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_beta_selic[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "da Selic"))
}

#par(old.par) 


#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_sd_selic[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "da Selic"))
}

#par(old.par) 

#####################################
# Para o IBC-Br
#####################################

#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_beta_ibcbr[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "do IBC-Br"))
}

#par(old.par) 


#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_sd_ibcbr[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "do IBC-Br"))
}

#par(old.par)

#####################################
# Para o Câmbio
#####################################

#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_beta_cambio[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "do Câmbio"))
}

#par(old.par) 


#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_sd_cambio[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "do Câmbio"))
}

#par(old.par)

#####################################
# Para o IPCA
#####################################

#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_beta_ipca[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "do IPCA"))
}

#par(old.par) 


#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
  cumuplot(mcmc_sd_ipca[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "do IPCA"))
}

#par(old.par)

Por fim, tem o diagnóstico do Geweke, que calcula uma estatística com distribuição normal e que compara o início e o final da distribuição. De novo vou plotar só uma parte para não ficar muito grande

geweke.plot(mcmc_beta_ct[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_beta_selic[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_beta_ibcbr[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_beta_cambio[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_beta_ipca[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_sd_ct[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_sd_selic[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_sd_ibcbr[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_sd_cambio[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

geweke.plot(mcmc_sd_ipca[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)

Diagnóstico da autocorrelação

O bloco abaixo eu montei para entender melhor o que a autocorr() do pacote coda faz. No fim, vou ficar com a função de autocorrelação (não a parcial) pois ela avalia o efeito até o instante j.

# Tentando descobrir o que as coisas são
# Esse é o do coda, ele é o mais parecido com a acf, mas não chega a ser igual.
autocorr(mcmc_sd_ct[,i], lags = c(20,100,1000))

c(acf(mcmc_sd_ct[,1], lag.max = 20, plot = F)[20]$acf,acf(mcmc_sd_ct[,1], lag.max = 100, plot = F)[100]$acf, acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F)[1000]$acf)
c(Acf(mcmc_sd_ct[,1], lag.max = 20, plot = F, type = "correlation")[20]$acf, Acf(mcmc_sd_ct[,1], lag.max = 100, plot = F, type = "correlation")[100]$acf, Acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F, type = "correlation")[1000]$acf)

# Esse não é parecido com ninguém
c(Acf(mcmc_sd_ct[,1], lag.max = 20, plot = F, type = "covariance")[20]$acf, Acf(mcmc_sd_ct[,1], lag.max = 100, plot = F, type = "covariance")[100]$acf, Acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F, type = "covariance")[1000]$acf)


# Esses 3 são iguais
# Function Pacf computes (and by default plots) an estimate of the partial autocorrelation function of a (possibly multivariate) time series
# Aviso sobre as funções Acf e Pacf: The functions improve the acf, pacf and ccf functions. The main differences are that Acf does not plot a spike at lag 0 when type=="correlation" (which is redundant) and the horizontal axes show lags in time units rather than seasonal units..
c(pacf(mcmc_sd_ct[,1], lag.max = 20, plot = F)[20]$acf,pacf(mcmc_sd_ct[,1], lag.max = 100, plot = F)[100]$acf, pacf(mcmc_sd_ct[,1], lag.max = 1000, plot = F)[1000]$acf)
c(Acf(mcmc_sd_ct[,1], lag.max = 20, plot = F, type = "partial")[20]$acf, Acf(mcmc_sd_ct[,1], lag.max = 100, plot = F, type = "partial")[100]$acf, Acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F, type = "partial")[1000]$acf)
c(Pacf(mcmc_sd_ct[,1], lag.max = 20, plot = F)[20]$acf, Pacf(mcmc_sd_ct[,1], lag.max = 100, plot = F)[100]$acf, Pacf(mcmc_sd_ct[,1], lag.max = 1000, plot = F)[1000]$acf)

Agora estou calculando a função te autocorrelação para depois fazer uns fancy graphs.

autocorrelacao_beta_ct <- vector()
autocorrelacao_beta_selic <- vector()
autocorrelacao_beta_ibcbr <- vector()
autocorrelacao_beta_cambio <- vector()
autocorrelacao_beta_ipca <- vector()

autocorrelacao_sd_ct <- vector()
autocorrelacao_sd_selic <- vector()
autocorrelacao_sd_ibcbr <- vector()
autocorrelacao_sd_cambio <- vector()
autocorrelacao_sd_ipca <- vector()

pm <- proc.time()
for (i in 1:dim(mcmc_beta_ct)[2]) {
    autocorrelacao_beta_ct <- cbind(autocorrelacao_beta_ct, autocorr(mcmc_beta_ct[,i], lags = c(20,100,1000)))
    autocorrelacao_beta_selic <- cbind(autocorrelacao_beta_selic, autocorr(mcmc_beta_selic[,i], lags = c(20,100,1000)))
    autocorrelacao_beta_ibcbr <- cbind(autocorrelacao_beta_ibcbr, autocorr(mcmc_beta_ibcbr[,i], lags = c(20,100,1000)))
    autocorrelacao_beta_cambio <- cbind(autocorrelacao_beta_cambio, autocorr(mcmc_beta_cambio[,i], lags = c(20,100,1000)))
    autocorrelacao_beta_ipca <- cbind(autocorrelacao_beta_ipca, autocorr(mcmc_sd_ipca[,i], lags = c(20,100,1000)))
    autocorrelacao_sd_ct <- cbind(autocorrelacao_sd_ct, autocorr(mcmc_sd_ct[,i], lags = c(20,100,1000)))
    autocorrelacao_sd_selic <- cbind(autocorrelacao_sd_selic, autocorr(mcmc_sd_selic[,i], lags = c(20,100,1000)))
    autocorrelacao_sd_ibcbr <- cbind(autocorrelacao_sd_ibcbr, autocorr(mcmc_sd_ibcbr[,i], lags = c(20,100,1000)))
    autocorrelacao_sd_cambio <- cbind(autocorrelacao_sd_cambio, autocorr(mcmc_sd_cambio[,i], lags = c(20,100,1000)))
    autocorrelacao_sd_ipca <- cbind(autocorrelacao_sd_ipca, autocorr(mcmc_sd_ipca[,i], lags = c(20,100,1000)))
}
proc.time() - pm

##    user  system elapsed 
##  150.36    0.25  152.46

Agora vamos para os gráficos…. a ideia é fazer algo como o Primiceri tem no anexo B (figura 9 - topo), porém dividindo por variável do VAR.

# Montando os data frames

#############################################
#############################################
# Betas
#############################################
#############################################
autocorrelacao_beta_ct <- data.frame(1:152, t(autocorrelacao_beta_ct))
names(autocorrelacao_beta_ct) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_ct <- melt(autocorrelacao_beta_ct, id.vars = "tempo")
autocorrelacao_beta_ct[,2] <- factor(autocorrelacao_beta_ct[,2])

autocorrelacao_beta_selic1 <- data.frame(1:152, t(autocorrelacao_beta_selic))
names(autocorrelacao_beta_selic1) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_selic1 <- melt(autocorrelacao_beta_selic1, id.vars = "tempo")
autocorrelacao_beta_selic1[,2] <- factor(autocorrelacao_beta_selic1[,2])

autocorrelacao_beta_ibcbr <- data.frame(1:152, t(autocorrelacao_beta_ibcbr)) 
names(autocorrelacao_beta_ibcbr) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_ibcbr <- melt(autocorrelacao_beta_ibcbr, id.vars = "tempo")
autocorrelacao_beta_ibcbr[,2] <- factor(autocorrelacao_beta_ibcbr[,2])


autocorrelacao_beta_cambio <- data.frame(1:152, t(autocorrelacao_beta_cambio))
names(autocorrelacao_beta_cambio) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_cambio <- melt(autocorrelacao_beta_cambio, id.vars = "tempo")
autocorrelacao_beta_cambio[,2] <- factor(autocorrelacao_beta_cambio[,2])

autocorrelacao_beta_ipca <- data.frame(1:152, t(autocorrelacao_beta_ipca))
names(autocorrelacao_beta_ipca) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_ipca <- melt(autocorrelacao_beta_ipca, id.vars = "tempo")
autocorrelacao_beta_ipca[,2] <- factor(autocorrelacao_beta_ipca[,2])

#############################################
#############################################
# Volatilidade
#############################################
#############################################
autocorrelacao_sd_ct <- data.frame(1:152, t(autocorrelacao_sd_ct))
names(autocorrelacao_sd_ct) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_ct <- melt(autocorrelacao_sd_ct, id.vars = "tempo")
autocorrelacao_sd_ct[,2] <- factor(autocorrelacao_sd_ct[,2])

autocorrelacao_sd_selic <- data.frame(1:152, t(autocorrelacao_sd_selic))
names(autocorrelacao_sd_selic) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_selic <- melt(autocorrelacao_sd_selic, id.vars = "tempo")
autocorrelacao_sd_selic[,2] <- factor(autocorrelacao_sd_selic[,2])

autocorrelacao_sd_ibcbr <- data.frame(1:152, t(autocorrelacao_sd_ibcbr)) 
names(autocorrelacao_sd_ibcbr) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_ibcbr <- melt(autocorrelacao_sd_ibcbr, id.vars = "tempo")
autocorrelacao_sd_ibcbr[,2] <- factor(autocorrelacao_sd_ibcbr[,2])

autocorrelacao_sd_cambio <- data.frame(1:152, t(autocorrelacao_sd_cambio))
names(autocorrelacao_sd_cambio) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_cambio <- melt(autocorrelacao_sd_cambio, id.vars = "tempo")
autocorrelacao_sd_cambio[,2] <- factor(autocorrelacao_sd_cambio[,2])

autocorrelacao_sd_ipca <- data.frame(1:152, t(autocorrelacao_sd_ipca))
names(autocorrelacao_sd_ipca) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_ipca <- melt(autocorrelacao_sd_ipca, id.vars = "tempo")
autocorrelacao_sd_ipca[,2] <- factor(autocorrelacao_sd_ipca[,2])

#############################################
#############################################
# Gráficos
#############################################
#############################################

##################### BETAS ########################

########################
# Razão capital trabalho
########################
p1a <- ggplot(autocorrelacao_beta_ct, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="Razão Cap./Trab.", lim = c(-0.1,0.6)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p1a <- p1a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", legend.text = element_text(size = 8), legend.title = element_text(size = 8), axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p1a <- p1a + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not use a gambiarra?
p1a <- p1a + labs(colour='Lags')

########################
# selic
########################
p2a <- ggplot(autocorrelacao_beta_selic1, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="Selic", lim = c(-0.1,0.6)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p2a <- p2a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p2a <- p2a + labs(linetype = "Lags") 
p2a <- p2a + labs(colour='Lags')

########################
# ibcbr
########################
p3a <- ggplot(autocorrelacao_beta_ibcbr, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="IBCBr", lim = c(-0.1,0.6)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p3a <- p3a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p3a <- p3a + labs(linetype = "Lags") 
p3a <- p3a + labs(colour='Lags')

########################
# cambio
########################
p4a <- ggplot(autocorrelacao_beta_cambio, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="Câmbio", lim = c(-0.1,0.6)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p4a <- p4a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p4a <- p4a + labs(linetype = "Lags") 
p4a <- p4a + labs(colour='Lags')

########################
# ipca
########################
p5a <- ggplot(autocorrelacao_beta_ipca, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="IPCA", lim = c(-0.1,0.6)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p5a <- p5a + theme(axis.text.x = element_text(angle=25, vjust = 2, size = 5, hjust = 1), axis.title.x = element_text(size = 6, vjust = 5), axis.title.y = element_text(size =6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.0001,0.1), "cm")) # cima,direita, baixo, esquerda)
p5a <- p5a + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not?
p5a <- p5a + labs(colour='Lags')

#grid.draw(rbind(ggplotGrob(p1a), ggplotGrob(p2a), ggplotGrob(p3a), ggplotGrob(p4a), ggplotGrob(p5a), size = "first"))

################## VOLATILIDADE ######################


########################
# Razão capital trabalho
########################
p1 <- ggplot(autocorrelacao_sd_ct, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=T, aes(linetype = variable))+
        scale_y_continuous(name="Razão Cap./Trab.", lim = c(-0.1,0.5)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", legend.text = element_text(size = 8), legend.title = element_text(size = 8), axis.text.y = element_text(size = 6), legend.key.size = unit(.5, "cm"), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm"))
p1 <- p1 + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not use a gambiarra?
p1 <- p1 + labs(colour='Lags')

########################
# selic
########################

p2 <- ggplot(autocorrelacao_sd_selic, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="Selic", lim = c(-0.1,0.5)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm"))
p2 <- p2 + labs(linetype = "Lags") 
p2 <- p2 + labs(colour='Lags')

########################
# ibcbr
########################

p3 <- ggplot(autocorrelacao_sd_ibcbr, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="IBCBr", lim = c(-0.1,0.5)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm"))
p3 <- p3 + labs(linetype = "Lags") 
p3 <- p3 + labs(colour='Lags')

########################
# cambio
########################

p4 <- ggplot(autocorrelacao_sd_cambio, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="Câmbio", lim = c(-0.1,0.5)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm")) # cima,direita, baixo, esquerda)
p4 <- p4 + labs(linetype = "Lags") 
p4 <- p4 + labs(colour='Lags')

########################
# ipca
########################

p5 <- ggplot(autocorrelacao_sd_ipca, aes(tempo, value, colour = variable)) +
        geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
        scale_y_continuous(name="IPCA", lim = c(-0.1,0.5)) +
        scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
        theme_bw() +
        scale_colour_brewer(palette="Set1")
        
p5 <- p5 + theme(axis.text.x = element_text(angle=25, vjust = 2, size = 5, hjust = 1), axis.title.x = element_text(size = 6, vjust = 5), axis.title.y = element_text(size =6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.0001,0.1), "cm")) # cima,direita, baixo, esquerda
p5 <- p5 + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not?
p5 <- p5 + labs(colour='Lags')



ggarrange(p1a, p1, p2a, p2, p3a, p3, p4a, p4, p5a, p5, common.legend = TRUE, ncol = 2, nrow = 5)

#################
# grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(p3), ggplotGrob(p4), ggplotGrob(p5), size = "first"))
# grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5) # This squish the graph with the legend

# Estimar outroS 4 modelos guardando apenas as últimas 200 observações e fazer uma anova não paramétrica

Diagnóstico multivariado

A estatística de Gelman-Rubin mede se há diferença significativa entre a variância dentro e entre diversas cadeias utilizando um fator de redução de escala. Aqui eu optei por simular duas cadeias e fazer a comparação entre elas. Em função da memória do computador eu utilizei um número menor de iterações salvas: ao invés de salvar todos os valores, armazenei apenas os 10.000 últimos valores (20%).

# Primeira coisa é limpar a memória
rm(list = ls())

# Auxiliary variables, so I don't need to bother when something changes
inicio <- "2003-02-01"
fim    <- "2017-10-31"

# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca   <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries

trabalho <- BETS.get("7620", from = inicio, to = fim)
capital  <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho

# Usando ndiffs(ibcbr,test="adf",alpha = 0.1) se encontra que o ibcbr é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
ibcbr_raw    <- BETS.get("24364", from = inicio_cambio, to = fim)
ibcbr        <- diff(log(ibcbr_raw), 1)

cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio     <- diff(log(cambio_raw), 1)

selic_4390 <- BETS.get("4390", from = inicio, to = fim) 
selic <- ((1+selic_4390/100)^(12)-1)*100

ipca_raw <- BETS.get("433", from = inicio_ipca, to = fim) 
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()

final <- length(ipca_acum)
for (i in 12:final){
  ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}

ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))

ipca <- ts(ipca,  start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD

m <- seas(x = capital_trabalho)
capital_trabalho2 <- final(m)

m <- seas(x = selic, transform.function = "none")
selic2 <- final(m)

var1 <- cbind(capital_trabalho2, selic2, ibcbr, cambio, ipca) # The bvar function does not allows data.frames

nburn. <- 5000
nrep. <- 50000
thinfac. <- 5

set.seed(6969)
pm <- proc.time()
fit1 <- bvar.sv.tvp(var1, p=1, tau = 24, nburn = nburn., nrep = nrep., nf = 1, thinfac = thinfac.)

## [1] "2018-01-26 10:01:39 -- now starting MCMC"
## [1] "2018-01-26 10:02:48 -- now at iteration 5000"
## [1] "2018-01-26 10:05:04 -- now at iteration 15000"
## [1] "2018-01-26 10:07:20 -- now at iteration 25000"
## [1] "2018-01-26 10:09:36 -- now at iteration 35000"
## [1] "2018-01-26 10:11:53 -- now at iteration 45000"
## [1] "2018-01-26 10:14:09 -- now at iteration 55000"

proc.time() - pm

##    user  system elapsed 
##  733.25   16.95  750.40

set.seed(9696)
fit2 <- bvar.sv.tvp(var1, p=1, tau = 24, nburn = nburn., nrep = nrep., nf = 1, thinfac = thinfac.)

## [1] "2018-01-26 10:14:10 -- now starting MCMC"
## [1] "2018-01-26 10:15:17 -- now at iteration 5000"
## [1] "2018-01-26 10:17:34 -- now at iteration 15000"
## [1] "2018-01-26 10:19:50 -- now at iteration 25000"
## [1] "2018-01-26 10:22:06 -- now at iteration 35000"
## [1] "2018-01-26 10:24:22 -- now at iteration 45000"
## [1] "2018-01-26 10:26:40 -- now at iteration 55000"

mcmc_1 <- mcmc(cbind(parameter.draws(fit1, type = "lag1", row = 1, col = 1), parameter.draws(fit1, type = "lag1", row = 2, col = 2), parameter.draws(fit1, type = "lag1", row = 3, col = 3), parameter.draws(fit1, type = "lag1", row = 4, col = 4), parameter.draws(fit1, type = "lag1", row = 5, col = 5)))

mcmc_2 <- mcmc(cbind(parameter.draws(fit2, type = "lag1", row = 1, col = 1), parameter.draws(fit2, type = "lag1", row = 2, col = 2), parameter.draws(fit2, type = "lag1", row = 3, col = 3), parameter.draws(fit2, type = "lag1", row = 4, col = 4), parameter.draws(fit2, type = "lag1", row = 5, col = 5)))

lista <- mcmc.list(mcmc_1, mcmc_2)
gelman.diag(lista)

## Potential scale reduction factors:
## 
##        Point est. Upper C.I.
##   [1,]       1.01       1.03
##   [2,]       1.01       1.03
##   [3,]       1.01       1.03
##   [4,]       1.01       1.03
##   [5,]       1.01       1.02
##   [6,]       1.01       1.03
##   [7,]       1.00       1.02
##   [8,]       1.00       1.02
##   [9,]       1.00       1.02
##  [10,]       1.00       1.02
##  [11,]       1.00       1.02
##  [12,]       1.00       1.01
##  [13,]       1.00       1.01
##  [14,]       1.00       1.01
##  [15,]       1.00       1.01
##  [16,]       1.00       1.01
##  [17,]       1.00       1.02
##  [18,]       1.00       1.01
##  [19,]       1.00       1.01
##  [20,]       1.00       1.01
##  [21,]       1.00       1.01
##  [22,]       1.00       1.01
##  [23,]       1.00       1.01
##  [24,]       1.00       1.02
##  [25,]       1.01       1.02
##  [26,]       1.01       1.03
##  [27,]       1.01       1.03
##  [28,]       1.01       1.03
##  [29,]       1.01       1.03
##  [30,]       1.00       1.02
##  [31,]       1.00       1.01
##  [32,]       1.00       1.01
##  [33,]       1.00       1.01
##  [34,]       1.00       1.01
##  [35,]       1.00       1.01
##  [36,]       1.00       1.01
##  [37,]       1.00       1.01
##  [38,]       1.00       1.02
##  [39,]       1.01       1.03
##  [40,]       1.01       1.03
##  [41,]       1.00       1.02
##  [42,]       1.00       1.01
##  [43,]       1.00       1.01
##  [44,]       1.00       1.01
##  [45,]       1.00       1.01
##  [46,]       1.00       1.00
##  [47,]       1.00       1.00
##  [48,]       1.00       1.00
##  [49,]       1.00       1.00
##  [50,]       1.00       1.00
##  [51,]       1.00       1.01
##  [52,]       1.00       1.01
##  [53,]       1.00       1.01
##  [54,]       1.00       1.00
##  [55,]       1.00       1.00
##  [56,]       1.00       1.00
##  [57,]       1.00       1.00
##  [58,]       1.00       1.00
##  [59,]       1.00       1.00
##  [60,]       1.00       1.01
##  [61,]       1.00       1.01
##  [62,]       1.00       1.00
##  [63,]       1.00       1.00
##  [64,]       1.00       1.00
##  [65,]       1.00       1.00
##  [66,]       1.00       1.01
##  [67,]       1.00       1.02
##  [68,]       1.01       1.04
##  [69,]       1.01       1.05
##  [70,]       1.01       1.04
##  [71,]       1.01       1.04
##  [72,]       1.01       1.04
##  [73,]       1.01       1.04
##  [74,]       1.01       1.04
##  [75,]       1.01       1.04
##  [76,]       1.01       1.04
##  [77,]       1.01       1.04
##  [78,]       1.01       1.03
##  [79,]       1.00       1.02
##  [80,]       1.00       1.00
##  [81,]       1.00       1.00
##  [82,]       1.00       1.00
##  [83,]       1.00       1.00
##  [84,]       1.00       1.00
##  [85,]       1.00       1.00
##  [86,]       1.00       1.00
##  [87,]       1.00       1.00
##  [88,]       1.00       1.00
##  [89,]       1.00       1.01
##  [90,]       1.00       1.00
##  [91,]       1.00       1.00
##  [92,]       1.00       1.00
##  [93,]       1.00       1.00
##  [94,]       1.00       1.00
##  [95,]       1.00       1.00
##  [96,]       1.00       1.00
##  [97,]       1.00       1.00
##  [98,]       1.00       1.00
##  [99,]       1.00       1.00
## [100,]       1.00       1.00
## [101,]       1.00       1.01
## [102,]       1.00       1.01
## [103,]       1.00       1.01
## [104,]       1.00       1.00
## [105,]       1.00       1.00
## [106,]       1.00       1.00
## [107,]       1.00       1.00
## [108,]       1.00       1.00
## [109,]       1.00       1.00
## [110,]       1.00       1.00
## [111,]       1.00       1.01
## [112,]       1.00       1.01
## [113,]       1.00       1.01
## [114,]       1.00       1.01
## [115,]       1.00       1.01
## [116,]       1.00       1.01
## [117,]       1.00       1.02
## [118,]       1.01       1.03
## [119,]       1.01       1.06
## [120,]       1.01       1.07
## [121,]       1.02       1.07
## [122,]       1.02       1.07
## [123,]       1.02       1.07
## [124,]       1.02       1.09
## [125,]       1.02       1.09
## [126,]       1.02       1.08
## [127,]       1.02       1.07
## [128,]       1.02       1.07
## [129,]       1.01       1.06
## [130,]       1.01       1.03
## [131,]       1.00       1.01
## [132,]       1.00       1.00
## [133,]       1.00       1.00
## [134,]       1.00       1.00
## [135,]       1.00       1.01
## [136,]       1.00       1.02
## [137,]       1.01       1.03
## [138,]       1.01       1.04
## [139,]       1.02       1.08
## [140,]       1.03       1.11
## [141,]       1.03       1.13
## [142,]       1.03       1.15
## [143,]       1.04       1.15
## [144,]       1.03       1.15
## [145,]       1.03       1.14
## [146,]       1.03       1.14
## [147,]       1.03       1.14
## [148,]       1.03       1.13
## [149,]       1.03       1.12
## [150,]       1.02       1.10
## [151,]       1.02       1.08
## [152,]       1.02       1.06
## [153,]       1.01       1.03
## [154,]       1.01       1.03
## [155,]       1.00       1.01
## [156,]       1.00       1.00
## [157,]       1.00       1.01
## [158,]       1.00       1.00
## [159,]       1.00       1.00
## [160,]       1.00       1.01
## [161,]       1.00       1.00
## [162,]       1.00       1.00
## [163,]       1.00       1.01
## [164,]       1.01       1.03
## [165,]       1.01       1.04
## [166,]       1.03       1.10
## [167,]       1.03       1.09
## [168,]       1.03       1.09
## [169,]       1.02       1.07
## [170,]       1.02       1.06
## [171,]       1.02       1.06
## [172,]       1.02       1.06
## [173,]       1.02       1.05
## [174,]       1.02       1.05
## [175,]       1.02       1.05
## [176,]       1.02       1.05
## [177,]       1.02       1.05
## [178,]       1.02       1.05
## [179,]       1.01       1.05
## [180,]       1.02       1.05
## [181,]       1.02       1.06
## [182,]       1.02       1.05
## [183,]       1.02       1.05
## [184,]       1.02       1.05
## [185,]       1.02       1.05
## [186,]       1.02       1.06
## [187,]       1.02       1.07
## [188,]       1.02       1.07
## [189,]       1.02       1.07
## [190,]       1.02       1.09
## [191,]       1.02       1.09
## [192,]       1.02       1.09
## [193,]       1.02       1.06
## [194,]       1.01       1.05
## [195,]       1.01       1.04
## [196,]       1.01       1.02
## [197,]       1.01       1.01
## [198,]       1.00       1.00
## [199,]       1.00       1.00
## [200,]       1.00       1.00
## [201,]       1.01       1.01
## [202,]       1.01       1.02
## [203,]       1.01       1.02
## [204,]       1.01       1.03
## [205,]       1.01       1.02
## [206,]       1.01       1.01
## [207,]       1.01       1.01
## [208,]       1.00       1.00
## [209,]       1.00       1.01
## [210,]       1.01       1.01
## [211,]       1.01       1.02
## [212,]       1.01       1.03
## [213,]       1.01       1.03
## [214,]       1.01       1.04
## [215,]       1.01       1.04
## [216,]       1.01       1.03
## [217,]       1.01       1.04
## [218,]       1.01       1.05
## [219,]       1.02       1.07
## [220,]       1.02       1.09
## [221,]       1.02       1.09
## [222,]       1.03       1.10
## [223,]       1.03       1.11
## [224,]       1.03       1.10
## [225,]       1.03       1.10
## [226,]       1.02       1.09
## [227,]       1.02       1.08
## [228,]       1.02       1.06
## [229,]       1.02       1.05
## [230,]       1.01       1.03
## [231,]       1.00       1.01
## [232,]       1.00       1.00
## [233,]       1.00       1.01
## [234,]       1.01       1.01
## [235,]       1.01       1.02
## [236,]       1.01       1.03
## [237,]       1.01       1.03
## [238,]       1.01       1.03
## [239,]       1.01       1.03
## [240,]       1.01       1.02
## [241,]       1.00       1.00
## [242,]       1.00       1.00
## [243,]       1.00       1.00
## [244,]       1.00       1.00
## [245,]       1.00       1.01
## [246,]       1.01       1.02
## [247,]       1.01       1.02
## [248,]       1.01       1.03
## [249,]       1.01       1.02
## [250,]       1.01       1.02
## [251,]       1.01       1.02
## [252,]       1.02       1.08
## [253,]       1.01       1.02
## [254,]       1.01       1.03
## [255,]       1.01       1.02
## [256,]       1.01       1.02
## [257,]       1.01       1.03
## [258,]       1.01       1.03
## [259,]       1.01       1.02
## [260,]       1.01       1.03
## [261,]       1.01       1.02
## [262,]       1.01       1.02
## [263,]       1.00       1.02
## [264,]       1.00       1.02
## [265,]       1.01       1.03
## [266,]       1.01       1.03
## [267,]       1.01       1.03
## [268,]       1.01       1.02
## [269,]       1.01       1.03
## [270,]       1.01       1.04
## [271,]       1.02       1.08
## [272,]       1.02       1.08
## [273,]       1.02       1.07
## [274,]       1.02       1.06
## [275,]       1.02       1.05
## [276,]       1.02       1.06
## [277,]       1.02       1.06
## [278,]       1.02       1.04
## [279,]       1.02       1.03
## [280,]       1.02       1.03
## [281,]       1.02       1.02
## [282,]       1.01       1.01
## [283,]       1.01       1.01
## [284,]       1.01       1.01
## [285,]       1.01       1.01
## [286,]       1.01       1.01
## [287,]       1.01       1.01
## [288,]       1.01       1.01
## [289,]       1.01       1.02
## [290,]       1.01       1.02
## [291,]       1.02       1.04
## [292,]       1.03       1.06
## [293,]       1.04       1.08
## [294,]       1.04       1.08
## [295,]       1.03       1.08
## [296,]       1.03       1.07
## [297,]       1.03       1.06
## [298,]       1.03       1.07
## [299,]       1.03       1.06
## [300,]       1.03       1.05
## [301,]       1.02       1.04
## [302,]       1.02       1.03
## [303,]       1.01       1.02
## [304,]       1.01       1.02
## [305,]       1.00       1.01
## [306,]       1.00       1.01
## [307,]       1.00       1.00
## [308,]       1.00       1.00
## [309,]       1.00       1.00
## [310,]       1.00       1.00
## [311,]       1.00       1.00
## [312,]       1.00       1.00
## [313,]       1.00       1.00
## [314,]       1.00       1.00
## [315,]       1.00       1.00
## [316,]       1.00       1.00
## [317,]       1.00       1.00
## [318,]       1.00       1.01
## [319,]       1.00       1.01
## [320,]       1.00       1.01
## [321,]       1.00       1.00
## [322,]       1.00       1.00
## [323,]       1.00       1.01
## [324,]       1.00       1.01
## [325,]       1.00       1.01
## [326,]       1.00       1.01
## [327,]       1.00       1.01
## [328,]       1.00       1.01
## [329,]       1.00       1.01
## [330,]       1.00       1.01
## [331,]       1.00       1.01
## [332,]       1.00       1.01
## [333,]       1.00       1.01
## [334,]       1.00       1.01
## [335,]       1.00       1.01
## [336,]       1.00       1.01
## [337,]       1.00       1.01
## [338,]       1.00       1.01
## [339,]       1.00       1.01
## [340,]       1.00       1.01
## [341,]       1.00       1.01
## [342,]       1.00       1.00
## [343,]       1.00       1.00
## [344,]       1.00       1.00
## [345,]       1.00       1.00
## [346,]       1.00       1.01
## [347,]       1.00       1.01
## [348,]       1.00       1.01
## [349,]       1.00       1.01
## [350,]       1.00       1.02
## [351,]       1.00       1.02
## [352,]       1.00       1.02
## [353,]       1.00       1.02
## [354,]       1.00       1.02
## [355,]       1.01       1.02
## [356,]       1.01       1.02
## [357,]       1.01       1.03
## [358,]       1.01       1.02
## [359,]       1.00       1.02
## [360,]       1.00       1.02
## [361,]       1.00       1.01
## [362,]       1.00       1.01
## [363,]       1.00       1.01
## [364,]       1.00       1.01
## [365,]       1.00       1.01
## [366,]       1.00       1.01
## [367,]       1.00       1.01
## [368,]       1.00       1.01
## [369,]       1.00       1.01
## [370,]       1.00       1.01
## [371,]       1.00       1.00
## [372,]       1.00       1.00
## [373,]       1.00       1.00
## [374,]       1.00       1.00
## [375,]       1.00       1.00
## [376,]       1.00       1.00
## [377,]       1.00       1.00
## [378,]       1.00       1.00
## [379,]       1.00       1.00
## [380,]       1.00       1.00
## [381,]       1.00       1.00
## [382,]       1.00       1.00
## [383,]       1.00       1.00
## [384,]       1.00       1.01
## [385,]       1.00       1.01
## [386,]       1.00       1.01
## [387,]       1.00       1.02
## [388,]       1.00       1.02
## [389,]       1.00       1.02
## [390,]       1.00       1.02
## [391,]       1.00       1.02
## [392,]       1.00       1.02
## [393,]       1.00       1.02
## [394,]       1.00       1.02
## [395,]       1.00       1.02
## [396,]       1.00       1.02
## [397,]       1.00       1.02
## [398,]       1.00       1.01
## [399,]       1.00       1.01
## [400,]       1.00       1.01
## [401,]       1.00       1.02
## [402,]       1.00       1.02
## [403,]       1.00       1.02
## [404,]       1.01       1.02
## [405,]       1.01       1.02
## [406,]       1.00       1.02
## [407,]       1.00       1.02
## [408,]       1.00       1.02
## [409,]       1.00       1.01
## [410,]       1.00       1.01
## [411,]       1.00       1.01
## [412,]       1.00       1.01
## [413,]       1.00       1.01
## [414,]       1.00       1.01
## [415,]       1.00       1.01
## [416,]       1.00       1.01
## [417,]       1.00       1.01
## [418,]       1.00       1.01
## [419,]       1.00       1.01
## [420,]       1.00       1.01
## [421,]       1.00       1.01
## [422,]       1.00       1.00
## [423,]       1.00       1.00
## [424,]       1.00       1.00
## [425,]       1.00       1.00
## [426,]       1.00       1.00
## [427,]       1.00       1.00
## [428,]       1.00       1.00
## [429,]       1.00       1.00
## [430,]       1.00       1.00
## [431,]       1.00       1.00
## [432,]       1.00       1.00
## [433,]       1.00       1.00
## [434,]       1.00       1.00
## [435,]       1.00       1.01
## [436,]       1.00       1.01
## [437,]       1.00       1.01
## [438,]       1.00       1.01
## [439,]       1.00       1.01
## [440,]       1.00       1.01
## [441,]       1.00       1.01
## [442,]       1.00       1.01
## [443,]       1.00       1.01
## [444,]       1.01       1.02
## [445,]       1.01       1.02
## [446,]       1.01       1.02
## [447,]       1.01       1.02
## [448,]       1.01       1.02
## [449,]       1.01       1.02
## [450,]       1.01       1.02
## [451,]       1.01       1.02
## [452,]       1.01       1.02
## [453,]       1.01       1.02
## [454,]       1.01       1.02
## [455,]       1.01       1.02
## [456,]       1.01       1.01
## [457,]       1.01       1.06
## [458,]       1.01       1.05
## [459,]       1.01       1.04
## [460,]       1.01       1.03
## [461,]       1.01       1.03
## [462,]       1.01       1.03
## [463,]       1.01       1.03
## [464,]       1.01       1.03
## [465,]       1.00       1.02
## [466,]       1.00       1.01
## [467,]       1.00       1.01
## [468,]       1.00       1.01
## [469,]       1.00       1.01
## [470,]       1.00       1.00
## [471,]       1.00       1.00
## [472,]       1.00       1.00
## [473,]       1.00       1.00
## [474,]       1.00       1.00
## [475,]       1.00       1.00
## [476,]       1.00       1.00
## [477,]       1.00       1.00
## [478,]       1.00       1.00
## [479,]       1.00       1.00
## [480,]       1.00       1.00
## [481,]       1.00       1.00
## [482,]       1.00       1.00
## [483,]       1.00       1.00
## [484,]       1.00       1.00
## [485,]       1.00       1.00
## [486,]       1.00       1.00
## [487,]       1.00       1.00
## [488,]       1.00       1.00
## [489,]       1.00       1.00
## [490,]       1.00       1.01
## [491,]       1.00       1.01
## [492,]       1.00       1.01
## [493,]       1.00       1.01
## [494,]       1.00       1.01
## [495,]       1.00       1.01
## [496,]       1.00       1.01
## [497,]       1.00       1.01
## [498,]       1.00       1.01
## [499,]       1.00       1.02
## [500,]       1.00       1.02
## [501,]       1.00       1.01
## [502,]       1.00       1.01
## [503,]       1.00       1.01
## [504,]       1.00       1.01
## [505,]       1.00       1.01
## [506,]       1.00       1.02
## [507,]       1.00       1.02
## [508,]       1.00       1.02
## [509,]       1.00       1.02
## [510,]       1.00       1.02
## [511,]       1.00       1.02
## [512,]       1.00       1.02
## [513,]       1.00       1.02
## [514,]       1.00       1.02
## [515,]       1.00       1.02
## [516,]       1.00       1.02
## [517,]       1.00       1.02
## [518,]       1.00       1.02
## [519,]       1.00       1.02
## [520,]       1.00       1.02
## [521,]       1.00       1.02
## [522,]       1.00       1.02
## [523,]       1.00       1.01
## [524,]       1.00       1.01
## [525,]       1.00       1.01
## [526,]       1.00       1.01
## [527,]       1.00       1.01
## [528,]       1.00       1.01
## [529,]       1.00       1.01
## [530,]       1.00       1.01
## [531,]       1.00       1.01
## [532,]       1.00       1.01
## [533,]       1.00       1.01
## [534,]       1.00       1.01
## [535,]       1.00       1.01
## [536,]       1.00       1.01
## [537,]       1.00       1.02
## [538,]       1.00       1.02
## [539,]       1.00       1.02
## [540,]       1.00       1.01
## [541,]       1.00       1.01
## [542,]       1.00       1.01
## [543,]       1.00       1.01
## [544,]       1.00       1.01
## [545,]       1.00       1.01
## [546,]       1.00       1.01
## [547,]       1.00       1.01
## [548,]       1.00       1.01
## [549,]       1.00       1.01
## [550,]       1.00       1.01
## [551,]       1.00       1.01
## [552,]       1.00       1.01
## [553,]       1.00       1.01
## [554,]       1.00       1.01
## [555,]       1.00       1.02
## [556,]       1.01       1.03
## [557,]       1.01       1.03
## [558,]       1.01       1.03
## [559,]       1.01       1.02
## [560,]       1.00       1.02
## [561,]       1.00       1.01
## [562,]       1.00       1.01
## [563,]       1.00       1.01
## [564,]       1.00       1.00
## [565,]       1.00       1.00
## [566,]       1.00       1.00
## [567,]       1.00       1.00
## [568,]       1.00       1.00
## [569,]       1.00       1.00
## [570,]       1.00       1.00
## [571,]       1.00       1.00
## [572,]       1.00       1.00
## [573,]       1.00       1.00
## [574,]       1.00       1.00
## [575,]       1.00       1.00
## [576,]       1.00       1.00
## [577,]       1.00       1.00
## [578,]       1.00       1.00
## [579,]       1.00       1.00
## [580,]       1.00       1.01
## [581,]       1.00       1.01
## [582,]       1.00       1.01
## [583,]       1.00       1.01
## [584,]       1.01       1.03
## [585,]       1.00       1.02
## [586,]       1.00       1.01
## [587,]       1.00       1.00
## [588,]       1.00       1.00
## [589,]       1.00       1.00
## [590,]       1.00       1.00
## [591,]       1.00       1.00
## [592,]       1.00       1.00
## [593,]       1.00       1.00
## [594,]       1.00       1.00
## [595,]       1.00       1.00
## [596,]       1.00       1.01
## [597,]       1.01       1.02
## [598,]       1.01       1.02
## [599,]       1.01       1.02
## [600,]       1.01       1.02
## [601,]       1.01       1.02
## [602,]       1.01       1.02
## [603,]       1.00       1.02
## [604,]       1.00       1.01
## [605,]       1.00       1.01
## [606,]       1.00       1.01
## [607,]       1.00       1.01
## [608,]       1.00       1.01
## [609,]       1.00       1.00
## [610,]       1.00       1.00
## [611,]       1.00       1.00
## [612,]       1.00       1.01
## [613,]       1.00       1.01
## [614,]       1.00       1.00
## [615,]       1.00       1.00
## [616,]       1.00       1.00
## [617,]       1.00       1.00
## [618,]       1.00       1.00
## [619,]       1.00       1.00
## [620,]       1.00       1.01
## [621,]       1.00       1.01
## [622,]       1.00       1.00
## [623,]       1.00       1.01
## [624,]       1.00       1.01
## [625,]       1.00       1.01
## [626,]       1.00       1.01
## [627,]       1.00       1.01
## [628,]       1.00       1.00
## [629,]       1.00       1.00
## [630,]       1.00       1.01
## [631,]       1.00       1.00
## [632,]       1.00       1.01
## [633,]       1.00       1.01
## [634,]       1.00       1.01
## [635,]       1.00       1.01
## [636,]       1.01       1.01
## [637,]       1.00       1.01
## [638,]       1.00       1.01
## [639,]       1.00       1.01
## [640,]       1.00       1.01
## [641,]       1.00       1.01
## [642,]       1.00       1.00
## [643,]       1.01       1.01
## [644,]       1.01       1.01
## [645,]       1.01       1.01
## [646,]       1.01       1.01
## [647,]       1.01       1.02
## [648,]       1.01       1.02
## [649,]       1.01       1.02
## [650,]       1.01       1.01
## [651,]       1.01       1.01
## [652,]       1.01       1.01
## [653,]       1.01       1.01
## [654,]       1.00       1.00
## [655,]       1.00       1.00
## [656,]       1.00       1.00
## [657,]       1.00       1.01
## [658,]       1.00       1.01
## [659,]       1.00       1.01
## [660,]       1.01       1.02
## [661,]       1.01       1.02
## [662,]       1.00       1.01
## [663,]       1.00       1.01
## [664,]       1.00       1.01
## [665,]       1.00       1.01
## [666,]       1.00       1.00
## [667,]       1.00       1.00
## [668,]       1.00       1.00
## [669,]       1.00       1.00
## [670,]       1.00       1.00
## [671,]       1.00       1.00
## [672,]       1.00       1.00
## [673,]       1.00       1.00
## [674,]       1.00       1.00
## [675,]       1.00       1.01
## [676,]       1.01       1.01
## [677,]       1.01       1.02
## [678,]       1.01       1.02
## [679,]       1.01       1.02
## [680,]       1.01       1.02
## [681,]       1.01       1.02
## [682,]       1.01       1.02
## [683,]       1.01       1.02
## [684,]       1.01       1.02
## [685,]       1.01       1.03
## [686,]       1.01       1.02
## [687,]       1.01       1.01
## [688,]       1.00       1.00
## [689,]       1.00       1.01
## [690,]       1.00       1.01
## [691,]       1.00       1.01
## [692,]       1.00       1.01
## [693,]       1.00       1.01
## [694,]       1.00       1.01
## [695,]       1.00       1.00
## [696,]       1.00       1.00
## [697,]       1.00       1.00
## [698,]       1.00       1.00
## [699,]       1.00       1.00
## [700,]       1.00       1.01
## [701,]       1.00       1.01
## [702,]       1.01       1.01
## [703,]       1.01       1.01
## [704,]       1.01       1.01
## [705,]       1.01       1.01
## [706,]       1.00       1.00
## [707,]       1.00       1.01
## [708,]       1.00       1.01
## [709,]       1.01       1.02
## [710,]       1.00       1.02
## [711,]       1.00       1.02
## [712,]       1.00       1.02
## [713,]       1.00       1.01
## [714,]       1.00       1.01
## [715,]       1.00       1.01
## [716,]       1.00       1.00
## [717,]       1.00       1.00
## [718,]       1.01       1.01
## [719,]       1.00       1.00
## [720,]       1.00       1.00
## [721,]       1.00       1.01
## [722,]       1.01       1.01
## [723,]       1.01       1.01
## [724,]       1.01       1.01
## [725,]       1.01       1.01
## [726,]       1.01       1.02
## [727,]       1.02       1.05
## [728,]       1.03       1.06
## [729,]       1.03       1.07
## [730,]       1.02       1.06
## [731,]       1.03       1.06
## [732,]       1.03       1.08
## [733,]       1.03       1.10
## [734,]       1.03       1.09
## [735,]       1.03       1.09
## [736,]       1.04       1.11
## [737,]       1.04       1.10
## [738,]       1.03       1.06
## [739,]       1.02       1.04
## [740,]       1.01       1.02
## [741,]       1.01       1.01
## [742,]       1.01       1.03
## [743,]       1.01       1.02
## [744,]       1.01       1.01
## [745,]       1.01       1.01
## [746,]       1.00       1.01
## [747,]       1.00       1.01
## [748,]       1.01       1.03
## [749,]       1.01       1.05
## [750,]       1.01       1.05
## [751,]       1.01       1.04
## [752,]       1.01       1.04
## [753,]       1.01       1.03
## [754,]       1.01       1.03
## [755,]       1.01       1.03
## [756,]       1.01       1.03
## [757,]       1.00       1.02
## [758,]       1.00       1.02
## [759,]       1.00       1.01
## [760,]       1.00       1.01
## 
## Multivariate psrf
## 
## 1.37

gelman.plot(lista)

# Tentar descobrir como colocar label no gelman plot